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Relativistic temperature of gas raises the issue of the equation of state (EoS) 
in relativistic hydrodynamics. We study the EoS for numerical relativistic hy- 
drodynamics, and propose a new EoS that is simple and yet approximates very 
closely the EoS of the single- component perfect gas in relativistic regime. We also 
discuss the calculation of primitive variables from conservative ones for the EoS's 
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■ 1. Introduction 

> ■ 

I Relativistic flows are involved in many high-energy astrophysical phenomena. Examples 

^ I includes relativistic jets from Galactic sources (see Mirabel Sz Rodriguez 1999, for reviews), 
extragalactic jets from AGNs (see Zensus 1997, for reviews), and gamma-ray bursts (see 
Meszaros 2002, for reviews). In relativistic jets from some Galactic microquasars, intrinsic 
beam velocities larger than 0.9c are typically required to explain the observed superluminal 
motions. In some powerful extragalactic radio sources, ejections from galactic nuclei produce 
true beam velocities of more than 0.98c. In the general fireball model of gamma-ray bursts, 
the internal energy of gas is converted into the bulk kinetic energy during expansion and this 
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expansion leads to relativistic outflows with high bulk Lorentz factors > 100. The flow mo- 
tions in these objects are usually highly nonlinear and intrinsically complex. Understanding 
such relativistic flows is important for correctly interpreting the observed phenomena, but 
often studying them is possible only through numerical simulations. 

Numerical codes for special relativistic hydrodynamics (hereafter RHDs) have been suc- 
cessfully built, based on explicit finite difference upwind schemes that were originally de- 
veloped for codes of non-relativistic hydrodynamics. These schemes utilize approximate or 
exact Riemann solvers and the characteristic decomposition of the hyperbolic system of 
conservation equations. RHD codes based on upwind schemes arc able to capture sharp 
discontinuities robustly in complex flows, and to describe the physical solution reliably. A 
partial hst of such codes includes the foUowings: Falle & Komissarov (1996) based on the 
van Leer scheme, Marti & Miiller (1996), Aloy et al. (1999), and Mignone et al. (2005) 
based on the PPM scheme, Sokolov et al. (2001) based on the Godunov scheme, Choi & Ryu 
(2005) based on the TVD scheme, Dolezal & Wong (1995), Donat et al (1998), DelZanna & 
Bucciantini (2002), and Rahman & Moore (2005) based on the ENO scheme, and Mignone 
& Bodo (2005) based on the HLL scheme. Reviews of some numerical approaches and test 
problems can be found in Marti & Miiller (2003) and Wilson & Mathews (2003). 

Gas in RHDs is characterized by relativistic fluid speed {v ~ c) and/or relativistic 
temperature (internal energy much greater than rest energy), and the latter brings us to the 
issue of the equation of state (hereafter EoS) of the gas. The EoS most commonly used in 
numerical RHDs, which is designed for the gas with constant ratio of specific heats, however, 
is essentially valid only for the gas of either subrelativistic or ultrarelativistic temperature. It 
is because that is not derived from relativistic kinetic theory. On the other hand, the EoS of 
the single-component perfect gas in relativistic regime can be derived from thermodynamics. 
But its form involves the modified Bessel functions (see Synge 1957), and is too comphcated 
to be implemented in numerical schemes. 

In this paper, we study EoS for numerical RHDs. We first revisit two EoS's previously 
used in numerical codes, specifically the one with constant ratio of specific heats, and the 
other first used by Mathews (1971) and later proposed for numerical RHDs by Mignone et al. 
(2005). Wc then propose a new EoS which is simple to be implemented in numerical codes 
with minimum efforts and minimum computational costs, but at the same time approximates 
very closely the EoS of the single- component perfect gas in relativistic regime. We also 
discuss the calculation of primitive variables from conservative ones for the three EoS's. 
Then we present the entire eigenstructure of RHDs for a general EoS, in a way to be used to 
build numerical codes. In order to see the consequence of different EoS's, shock tube tests 
performed with a code based on the TVD scheme are presented. The tests demonstrate 



-3- 



the differences in fiow structure due to different EoS's. Employing a correct EoS sliould 
be important to get quantitatively correct results in problems involving a transition from 
non-relativistic temperature to relativistic temperature or vice versa. 

This paper is organized as follows. In sections 2 and 3 we discuss three EoS's and the 
calculation of primitive variables from conservative ones for those three. In sections 4 we 
present the cigcnstructure of RHDs with a general EoS. In sections 5 and 6 we present a 
code based on the TVD scheme and shock tube tests with the code. Concluding remarks 
are drawn in section 7. 



2. Relativistic Hydrodynamics 

2.1. Bcisic Equations 

The special RHD equations for an ideal fluid can be written in the laboratory frame of 
reference as a hyperbolic system of conservation equations 

^ + ^W"^+I»5«) = 0, (16) 

^ + ^[(£ + p)«,l=0. (ic) 

where Mj, and E arc the mass density, momentum density, and total energy density, 
respectively (see, e.g., Landau & Lifshitz 1959; Wilson & Mathews 2003). The conserved 
quantities in the laboratory frame are expressed as 

D = rp, (2a) 

Mi = Vphvi, (26) 

E = T^ph - p, (2c) 

where p, Vi, p, and h are the proper mass density, fluid three- velocity, isotropic gas pressure 
and specific enthalpy, respectively, and the Lorentz factor is given by 

r^^=L= with v^^vl + vl + vl (3) 



V 



In the above, the Latin indices {e.g., i) represents spatial coordinates and conventional 
Einstein summation is used. The speed of light is set to unity (c = 1) throughout this paper. 



-4- 



2.2. Equation of State 

The above system of equations is closed with an EoS. Without loss of generality it is 
given as 

h=h{p,p). (4) 

Then the general form of polytropic index, n, and the general form of sound speed, Cg, 
respectively can be written as 

_ dh ^ ^ _ p dh , , 

^ ^ dp ' nh dp 

In addition we use a variable 7/1 to present the EoS property conveniently, 

Ih = (6) 

where = p/p is a tempcrature-like variable. 

The most commonly used EoS, which is called the ideal EoS (hereafter ID), is given as 

p-(7-l)(e-p) or h = l + -^ (7) 

7-1 

with a constant 7. Here 7 = Cp/ Cy is the ratio of specific heats, and e is the sum of the internal 
and rest-mass energy densities in the local frame and is related to the specific enthalpy as 

..t±l. ,S) 

For ID, 7/j = 7/(7 — 1) docs not depend on 0. ID may be correctly applied to the gas of 
either subrelativistic temperature with 7 = 5/3 or ultrarelativistic temperature with 7 = 4/3. 
But ID is rented from non-relativistic thermodynamics, and hence it is not consistent with 
relativistic kinetic theory. For example, we have 

^ 7-1' 70 + 7-1- 

In the high temperature limit, i.e., 0— >oo, and for 7 > 2, Cg > 1 i.e., admits superluminal 
sound speed. More importantly, using relativistic kinetic theory Taub (1948) showed that 
the choice of EoS is not arbitrary and has to satisfy the inequality, 

(/?,-0)(/i-40) > 1. (10) 



This rules out ID for 7 > 4/3, if applied for < < 00. 
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The correct EoS for the single-component perfect gas in relativistic regime (hereafter 
RP) can be derived (see Synge 1957), and is given as 

where K2 and are the modified Bessel functions of the second kind of order two and 
three, respectively. In the non- relativistic temperature limit (© 0), 7/j — > 5/2, and 
in the ultrarelativistic temperature limit (© — > 00), 4. However, using the above 

EoS comes with a price of extra computational costs (Falle & Komissarov 1996), since the 
thermodynamics of the fluid is expressed in terms of the modified Bessel functions. 

There have been efforts to find approximate EoS's which are simpler than RP but more 
accurate than ID. For example, Sokolov et al. (2001) proposed 



e = - (^/i - - j or /i = 26 + ^492 + 1. (12) 

But this EoS does not satisfy either Taub's inequality nor is consistent with the value of 7/1 
in the non- relativistic temperature limit. 

In a recent paper, Mignone et al. (2005) proposed for numerical RHDs an EoS that fits 
RP well. The EoS, which was first used by Mathews (1971), is given as 



1 = l(^__p\ or /. = ^e + ^Je2 + l, (13) 

p 3 Vp e) 2 2V 9' ^ ' 

and is abbreviated as TM following Mignone et al. (2005). With TM the expressions of n 
and Cs become 

, 3 e 5ev/e^T479 + 3e2 

2 2^92 + 4/9' ' 120792 + 4/9 + 1202 + 2 

TM corresponds to the lower bound of Taub's inequality, i.e., {h — Q){h — 40) = 1. It 
produces the right asymptotic values for jh- 

In this paper we propose a new EoS, which is a simpler algebraic function of and is 
also a better fit of RP compared to TM. We abbreviate our proposed EoS as RC and give it 
by 

p 3p + 2p 60^ + 40 + 1 

= or h = 2 . 15 

e-p 9p + 3p 30 + 2 ^ ^ 

With RC the expressions of n and Cg become 

_ 90^ + 120 + 2 2 _ e(30 + 2) (180^ + 240 + 5) 
^ ~ (30 + 2)2 ' ~ 3(602 + 40 + 1)(902 + 120 + 2) ' ^ ^ 
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RC satisfies Taub's inequality, {h — Q){h — 40) > 1, for all ©. It also produces the right 
asymptotic values for 7/j. For both TM and RC, we have correctly — > 50/3 in the 
non-relativistic temperature limit and 1/3 in the ultrarelativistic temperature limit, 

respectively. 

In Figure 1, ^h, and Cg are plotted as a function of © to compare TM and RC to RP 
as well as ID. One can see that RC is a better fit of RP than TM with 

\hTM - hnp\ ^ \hKC - hnp\ < Q_g^^_ ^^^^ 



^RP ^ ^RP 

It is to be remembered that both 7^ and n are independent of ©, if ID is used 



3. Calculation of Primitive Variables 

The RHD equations evolve the conserved quantities, D, Mj and E, but we need to 
know the values of the primitive variables, p, Vi, p, to solve the equations numerically. The 
primitive variables can be calculated by inverting the equations (2a-2c). The equations (2a- 
2c) explicitly include h, and here we discuss the inversion for the EoS's discussed in section 
2.2, that is, ID, TM, and RC. 



3.1. ID 

Schneider et al. (1993) showed that the equations (2a-2c) with the EoS in (7) reduce to 
a single quartic equation for v 



where 



+ biv^ + h2v' + 63^^ + 64 = 0, (18) 



_ 27(7 -l)Mi? _ 7^£;^ + 2(7-l)M^-(7-l)^/^^ , . 

' (7- 1)2(M2 + L>2)' ''2 (7 - 1)2(M2 + L>2) ' ^^"^ 



2 



(7 - 1)2(M2 + L»2) ' - (7 - 1)2(M2 + L»2) ' (19^) 

and M — + + M|. The quartic equation (18) can be solved numerically or analyt- 

ically. In Choi & Ryu (2005) the analytical solution was used for the very first time, though 
the exact nature of the solution was not presented. 

The general form of analytical roots for quartic equations can be found in Abramowitz & 
Stegun (1972) or on webs such as "http://mathworld.wolfram.com/QuarticEquation.html". 
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One may even use softwares such as Mathematica or Maxima to find the roots. We found 
that out of the four roots of the quartic equation (18), two are complex and two are real. 
The two real roots are 

-B + - AC -B - VB^ - AC 

where 

B^-{h + y/bj-Ab2 + Ax,), C^-{x,-^xl-Ah), (21a) 
x^^{R + T"^)-s + {R-T-2ys-^, (216) 

= ^^i^^i^lp^, S= ^"'^"\ T = R^ + S\ (21c) 

oi = -62, 0,2 = 6163 - 464, 03 = 46264 - ^'3 - b\hA. {21d) 

Among the two real roots, the first one is the solution that satisfies the upper and lower 
limits imposed by Schneider et al. (1993), thus v — z^. Once v is found, the quantities p, Vi, 
p, are calculated by 

(22a) 

My M, 
"^x — -nrv, Vy — -—T-v, Vz — -ttV, (22b) 

p = (7 - 1) [{E - M,v, - MyVy - M,v,) - p\ . (22c) 



3.2. TM 

Combining the equations (2a-2c) with the EoS in (13), we get a cubic equation for 
W = T^ -I 

+ c^W'^ + C2W + C3 = 0, (23) 

where 

_ {E'' + M^)[4(£;^ + W) - + £>^)] - lAWE'' 

~ 2{E^ - M2)2 ' ^2^"^ 

_ [4(£;2 + M^) - (M^ + L>2)]2 _ 57m2£;2 
- 16(£;2 - M-^Y ' ^^^^^ 

9m2£;2 

^^ = -16(£;^-M^)^- ^24c) 
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Cubic equations admit analytical solutions simpler than quartic equations (see also 
Abramowitz & Stegun 1972). We found that out of the three roots of the cubic equation 
(23), two are unphysical giving F < 1, and only one gives the physical solution, which is 

W = 2v^cos (0 - |, (25) 

where 

, 3c2-c2 H 9ciC2 - 27c3 - 2c? 

J= — -, cost = —z , if = — -. (25) 

9 ' 54 ^ ^ 

Then the fluid speed is calculated by 

W 

V = , (27) 
and the quantities p, Vi, p, are calculated by 

p=-. (28a) 

^oQM 

M ' ^ M ' M ' ^ ^ 

^ {E - M,v, - MyVy - M,v,f - 



3.3. R,C 

Combining the equations (2a-2c) with the EoS in (15), we get 

MVr2 - 1 [3Er{sr^ - 1) + 2r>(i - 4r2)] 

= 3T^ [4(M2 + £;2)r2 - (M^ + 4£;2)] - 2D{AET - D){T^ - 1). (29) 
Further simplification reduces it into an equation of 8*^ power in F. 

Although the equation (29) has to be solved numerically, it behaves very well. We first 
analyzed the nature of the roots with a root-finding routine in the IMSL library. As noted 
by Schneider et al. (1993), the physically meaningful solution should be between the upper 
hmit, F„, 

1 M 

V 1 - ^ 

and the lower limit, F;, that is derived inserting D = into equation (29): 
16(M2 - E^frf - 8(M2 - E^){M^ - AE^)rf + (M^ - 9M^E^ + 16E'^)r^ + M^E^ = (31) 
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(a cubic equation of Tf). Out of the eight roots of the equation (29), four are complex and 
four are real. Out of the four real roots, two are negative and two are positive. And out of 
the two real and positive roots, one is always larger than F^, and the other is between Fi 
and and so is the physical solution. 

Inside RHD codes the physical solution of equation (29) can be easily calculated by the 
Ncwton-Raphson method. With an initial guess F = or any value smaller than it including 
1, iteration can be proceeded upwards. Since the equation is extremely well-behaved, the 
iteration converges within a few steps. Once F is known, the fluid speed is calculated by 



and the quantities p, Vi, p, are calculated by 



V = ^^-^^ (32) 



Mx My M, 

Vx — inrV, Vy — —rV, — irrV (66b) 

_ {E - MjVi) -2p+ [{E - MjVif + 4p{E - M^v,) - 4p^]^ 

P — g ) {66c) 

where 

M,v^ = M^Vx + MyVy + M,v,. (34) 



4. Eigenvalues and Eigenvectors 

In building a code based on the Roe-type schemes such as the TVD and ENO schemes 
that solves a hyperbolic system of conservation equations, the eigenstructure (eigenvalues 
and eigenvectors of the Jacobian matrix) is required. The Eigenstructure for RHDs was 
previously described, for instance, in Donat et al. (1998). However, with the parameter 
vector different from that of Donat et al. (1998), the eigenvectors become different. Here we 
present our complete set of eigenvalues and eigenvectors without assuming any particular 
form of EoS. 



Equations (la)-(lc) can be written as 



dt dxj 



^^'■=0 (35) 
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with the state and flux vectors 





" D ' 




Dvj 


— * 






MiVj + pSij 




. E . 




. {E + p)vj _ 



(36) 



or as 



dt ~^ dxn 



0, 



^3 = 



(37) 



Here Aj is the 5x5 Jacobian matrix composed with the state and flux vectors. The con- 
struction of the matrix Aj can be simplified by introducing a parameter vector, u, as 



^ ^ dFj du 
^ du dq 

We choose the vector made of primitive variables as the parameter vector 



(38) 



u 



P 

Vi 

IP 



(39) 



4.1. One Velocity Component 



The eigenstructure is simplified if only a single component of velocity is chosen, i.e., 
V = Vx. In principle it can be reduced from that with three components of velocity in the 
next subsection. Nevertheless we present it, for the case that the simpler eigenstructure with 
one velocity component can be used. 



The explicit form of the Jacobian matrix. A, is presented in Appendix A. The eigenvalues 
of A are, 

V — Cs 

a- = , Co = V, 



The right eigenvectors are 



1 

Vh{v - c) 

.r/i(i - Csv)_ 



1 - CM 



Rq 



a+ 



1 + CgV 



Thv{l - nc2) 
. r/i(l - ncl) 



Th{v + c) 

rh(i + csv) 



and the left eigenvectors are 

1 



2hncl 



[h{l-nc1), r{v + ncs), - r{l + nCgv)] , 



(40) 



(41) 



(42a) 
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^ [h, rv, -r], (426) 



hnc^ 

^+--7Tr-^[H^-ncl), r{v-ncs), - r(l - nc,^;)] . (42c) 



2/inCg 

Here n and are given in equation (5). 



4.2. Three Velocity Components 



The a;-component of the Jacobian matrix, A^, when all the three components of velocity 
are considered, is presented in Appendix B. The eigenvalues of are 



{i-c^,)v,-Cs/r-VQ 



1 — c'^v'^ 



02 = Vx: 

03 = Vx, 



Os 



04 = V^, 



(43a) 

(436) 
(43c) 
(43ci) 

(43e) 



where Q 



cj.{Vy + v1). The eigenvalues represent the five characteristic speeds 



associated with two sound wave modes (ai and 05) and three entropy modes (02, 03, and 
04). A remarkable feature is that the eigenvalues do not explicitly depend on h and n, but 
only on Vi and Cg. Hence the eigenvalues are the same regardless of the choice of EoS once 
the sound speed is defined properly. 

The corresponding right eigenvectors {A^^R = aR), however, depends explicitly on h 
and n, and the complete set is given by 



Ri 



R2 — X [Xi, X2, X3, X4, ^^5]"^ , 



i?3 



Ra 



l-vl 



l-vl 



^, 2Vo,Vy, 1-vl + Vy, VyV^, 2Vy 



V. 

fh 



Rh 



1 - a5Vx 



, a5h{l-vl), h{l - a5Vx)vy, h{l - a5Vx)vz, h{l - vl) 



(44a) 

(446) 
(44c) 

{Ud) 
(44e) 
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where 



nci{v^y + vl) + {l-vl) 
Th 



X, = [2ncl{vl + vl) + (1 - ncl){l - vD] v 
Xs^ [nclivl + vl) + {l-vl)]vy, 
[ncliv^y + vl) + {1 - vD] V,, 
X, = 2ncl{vl + vl) + (1 - ncl){l - v^). 



X 



ncjil-vD' 



(45a) 

(456) 
(45c) 
(45d) 

(45e) 

(45/) 



The complete set of the left eigenvectors {LA^ — aL), which are orthonormal to the 
right eigenvectors, is 



Y 



[^11, Yl2, Yi3, Yi3, ¥15] , 



Lo = 



h 

p) "^xj '^y-i "^zj 



where 



L3 = [-Thvy, 0, 1, 0, 0] , 

L^^[-rhv„ 0, 0, 1, 0], 
1 

L5 — [Y51, Y52, ^53, Y53, Y55] , 

y 5 

h 

Yii = -p(l - aiVx){'^ - ncl), 
Yi2 = nai{l - clv"^) + 0,(1 + ncl)vl - (1 + n)vx 

YiS = -(1 + ^Cs)(l - flf^x-)^?;, 

Yi4 = -(1 + nc3)(l - aiV^)V;„ 
Yi5 = (1 + nc^v^) + (1 - c^)nv^ - aj(l + n)v^. 



K = /in 



(46a) 

(46&) 

(46c) 
(46d) 

(46e) 

(47a) 

(476) 

(47c) 
{47d) 
(47e) 

(47/) 



and index i — 1, 5. 



We note that with three degenerate modes that have same eigenvalues, a2 = as = a4, we 
have a freedom to write down the right and left eigenvectors in a variety of different forms. 
We chose to present the ones that produce the best results with the TVD code described 
next. 
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5. One-Dimensional Functioning Code 



To be used for demonstration of the differences in flow structure due to different EoS's, 
a one-dimensional functioning code based on the Total Variation Diminishing (TVD) scheme 
was built. The code utilizes the eigenvalues and eigenvectors given in the previous section, 
and can employ arbitrary EoS's including those in section 2.2. 



The TVD scheme, originally developed by Harten (f 983), is an Eulerian, finite-difference 
scheme with second-order accuracy in space and time. The second-order accuracy in time is 
achieved by modifying numerical flux using the quantities in five grid cells (see below and 
Harten 1983, for details). The scheme is basically identical to that previously used in Ryu et 
al. (1993) and Choi & Ryu (2005). But for completeness, the procedure is concisely shown 
here. 

The state vector at the cell center i at the time step n is updated by calculating the 
modified flux vector f x,i±i/2 along the x-direction at the cell interface i ± 1/2 as follows: 



5.1. 



The TVD Scheme 




(48) 




(49) 



fe=i 




At 



'fe,i+i/2 + ik,i+i/2)oik,i+i/2 - {gk,i + gk,i+l) , 



(50) 



Ax 



a 




(52a) 



(51) 



2|^fe,i+i/2|,2sign(^fc,i+i/2)^fe,j-i/2]}, 



(526) 



gk,i = sign(5(fc,i+i/2)max{0, mm[\gk,i+i/2\,2sign{gk,i+i/2)gk,i-i/2], 



min[2|5(fe,i+i/2|,sign(^fe,i+i/2)^fe,i_i/2]}, 



(52c) 




(53) 
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«fc,^+l/2 = £fe,^+l/2 • (?r+l - , (54) 

/ xV4efe + £fe for |a;| < 2£fe, 

L fI for \x\ > 2ek. 

Here, /c = 1 to 5 stand for the five cfiaracteristic modes. Tfie internal parameters e^s 
implicitly control numerical viscosity, and are defined for < Ek < 0.5. The flux limiters 
in equations (52a)-(52c) are the min-mod, monotonized central difference, and superbee 
limiters, respectively, a partial list of the limiters that are consistent with the TVD scheme, 
and one of them has to be employed. 



5.2. Quantities at Cell Interfaces 

To calculate the fluxes we need to define the local quantities at the cell interfaces, 
i + 1/2. The TVD scheme originally used the Roe's linearizion technique (Roe 1981) for 
it. Although it is possible to implement this linearizion technique in the relativistic domain 
in a computationally feasible way (see Eulderink & Mellema 1995), there is unlikely to be 
a significant advantage from the computational point of view. Instead, we simply use the 
algebraic averages of quantities at two adjacent cell centers to define the fluid three- velocity 
and specific enthalpy at the cell interfaces: 

_ + Vx,i+1 _ Vy,i + Vy,i+l _ Vz,i + Vz,i+1 /^-^n 
Vx,i+l/2 — ^ , — ^ , Vz,i+l/2 — ^ , [OO) 

hi + hi+i , , 

hi+i/2 = ^ ■ (57) 

Defining n and for the calculation of eigenvalues and eigenvectors at the cell interfaces 
depends on EoS. For ID, n is constant and 

For TM, we first compute from equation (13) 



5/ii+i/2 - ^9/i2^i/2 + 16 



e>iH-i/2 = '-^ , (59) 

then define ni+1/2 and Cs,i+i/2 according to equation (14). For RC, we first compute from 
equation (15) 



3hi+y2 - 8 + . /9/i2 ^ + 48/ii+i/2 - 32 
= (60) 



then define nj+1/2 and 05^^+1/2 according to equation (16). 
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6. Numerical Tests 

The differences induced by different EoS's are illustrated through a series of shock tube 
tests performed with the code described in the previous section. We use the tests used in 
previous works {e.g., Marti & Miiller 2003; Mignone et al. 2005), instead of inventing our own. 
Two sets are considered, one being purely one-dimensional with only the velocity component 
parallel to structure propagation, and the other with transverse velocity component. 

For the first set with parallel velocity component only, two tests are presented: 
PI: PL = 10, PR = 1, PL = 13.3, PR = 10~^, and Vp^L = Vp,R = initially, and tend = 0.45, 
P2: pl^ pr^ 1, PL = 10^, PR = 10~^, and Vp^L = Vp^R = initially, and tend = 0.4. 
The box covers the region of < x < 1. Here the subscripts L and R denote the quantities 
in the left and right states of the initial discontinuity at x — 0.5, and tend is the time when 
the solutions are presented. These two tests have been extensively used for tests of RHD 
codes with the ID EoS (see Marti & Miiller 2003), and the analytic solutions were described 
in Marti & Miiller (1994). 

Figures 2 and 3 show the numerical solutions with RC and TM, and the analytic solu- 
tions with ID and 7 = 5/3 and 4/3. The numerical solutions with RC and TM were obtained 
using the version of the TVD code having one velocity component (see section 4.1), and the 
analytic solutions with ID comes from the routine described in Marti & Miiller (1994). The 
numerical solutions with ID are almost indistinguishable from the analytic solutions, once 
they are calculated. 

The ID solutions with 7 = 4/3 and 5/3 show noticeable differences. The density shell 
between the contact discontinuity (hereafter CD) and the shock becomes thinner and taller 
with smaller 7, because the post shock pressure is lower and so is the shock propagation 
speed. The rarefaction wave is less elongated with 7 = 4/3, because the sound speed is 
lower. Those solutions with ID are also clearly different from the solutions obtained with 
RC and TM. The ID solution with 7 = 4/3 better approximates the solutions with RC and 
TM in the left region of the CD, because the flow has relativistic temperature of © > 1 there. 
The difference is, however, obvious in the shell between the CD and the shock, because G ~ 1 
there. On the other hand, the solutions obtained with RC and TM look very much alike. 
It reffects the similarity in the distributions of specific enthalpy in equations (13) and (15). 
Yet there is a noticeable difference, especially in the shell between the CD and the shock, 
and the difference in density reaches up to ~ 5%. 

For the second set with transverse velocity component, four tests, where different trans- 
verse velocities were added to the test P2, are presented: 
Tl: initially Vf^R = 0.99 to the right state, tend = 0.45, 
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T2: initially Vt,L — 0.9 to the left state, tend = 0.55, 

T3: initially Vt^L = Vt^R = 0.99 to the left and right states, tend = 0.18, 

T4: initially Vt^L = 0.9 and Vt^R = 0.99 to the left and right states, tend = 0.75. 

The notations are the same ones used in PI and P2. These are subsets of the tests originally 

suggested by Pons et al. (2000) with the ID EoS and later used by Mignone et al. (2005). 

Figures 4, 5, 6 and 7 show the numerical solutions with RC and TM and the analytic 
solutions with ID and 7 = 5/3 and 4/3. The numerical solutions with RC and TM were 
obtained using the version of the TVD code having three velocity components (see section 
4.2), and the analytic solutions with ID comes from the routine described in Pons et al. 
(2000). 

Again the ID solutions with 7 = 4/3 and 5/3 show noticeable differences. Especially 
with transverse velocity initially on the left side of the initial discontinuity (Figure 5, 6 and 
7), the parallel velocity reaches lower values, while the transverse velocity achieves higher 
values, with higher 7 = 5/3 in the region to the left of the CD. As a result, the density shell 
between the CD and the shock has propagated less. As in the P tests, the solutions with ID 
are clearly different from the solutions obtained with RC and TM, most noticeably in the 
shell between the CD and the shock. The solutions with RC and TM look very much ahke 
with differences in the density in the shell between the CD and the shock of about ~ 5%. 

We note that this paper is intended to focus on the EoS in numerical RHDs, not intended 
to present the performance of the code. Hence, one-dimensional tests of high resolution (with 
2^^ grid cells for the P tests and 2^^ grid cells the T tests) are presented to manifest the 
difference induced by different EoS's. The performance of the code such as capturing of 
shocks and CDs will be discussed elsewhere. 



7. Summary and Discussion 

The conservation equations for both Newtonian hydrodynamics and RHDs are strictly 
hyperbolic, rendering the apt use of upwind schemes for numerical codes. The actual imple- 
mentation to RHDs is, however, complicated, partly due to EoS. In this paper we study three 
EoS's for numerical RHDs, two being previously used and the other being newly proposed. 
The new EoS is simple and yet approximates the enthalpy of single-component perfect gas 
in relativistic regime with accuracy better than 0.8%. Then we discuss the calculation of 
primitive variables from conservative ones for the EoS's considered. We also present the 
eigenvalues and eigenvectors of RHDs for a general EoS, in a way that they are ready to be 
used to build numerical codes based on the Roe-type schemes such as the TVD and ENO 
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schemes. Finally we present numerical tests to show the differences in flow structure due to 
different EoS's 

The most commonly used, ideal EoS, can be used for the gas of entirely non-relativistic 
temperature (© <^ 1) with 7 = 5/3 or for the gas of entirely ultrarelativistic temperature 
(© ^ 1) with 7 = 4/3. However, if the transition from non-relativistic to relativistic or 
vice versa with ~ 0.1 — 1 is involved, the ideal EoS produces incorrect results and its use 
should be avoided. The EoS proposed by Mignone et al. (2005), TM, produces reasonably 
correct results with error of a few percent at most. The most preferable advantage of using 
TM is that the calculation of primitive variables admits analytic solutions, thereby making 
its implementation easy. The newly suggested EoS, RC, which approximates the EoS of the 
relativistic perfect gas, RP, most accurately, produces thermodynamically the most accurate 
results. At the same time it is simple enough to be implemented to numerical codes with 
minimum efforts and minimum computational costs. With RC the primitive variables should 
be calculated numerically by an iteration method such as the Newton-Raphson method. 
However, the equation for the calculation of primitive variables behaves extremely well, so 
the iteration converges in a few step without any trouble. 

In Galactic and extragalactic jets and gamma-ray bursts, as the flows travel relativistic 
fluid speeds (v ~ 1 but -C 1), they would hit the surrounding media. Then shocks are 
produced and the gas can be heated up to > 1. These kind of transitions, continuous or 
discontinuous, between relativistic bulk speeds and relativistic temperatures are intrinsic in 
astrophysical relativistic flows, and so a correct EoS is required to simulate the flows correctly. 
The correctness as well as the simplicity make RC suitable for astrophysical applications like 
these. 
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^2i = -p^(l-nc^) 
A22 = -p (1 - nc^) + 2vhn{l - c^) 
= -v^hn{l - c^) + A 

A32 = hn{l - c^v'^) 
N = hn(l - clv^) 



B. Jacobian Matrix with Three Velocity Components 
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^25 = -vlh{l + n) + h{l + nclv'^) 

^31 = - ncl) 

A32 = Vyh[n(l - cy) + vl{l + ncl)] 
A33 = v^h[n(l - cy) +vl{l + ncl)] 
A34 = v^VyV^h{l + ncl) 
Mb = -v^Vyh{l + n) 

Ml = ^V:,vji?'{\ - ncl) 
A42 = v,h[n{l - cy) + vl{l + ncl)] 

Ms = VxVyVzh{l + ncl) 
Au = v^h[n{l - cy) + vl{l + ncl)] 
^45 = -V:cVzh{l + n) 
A52 = hn{l - cy) 
N = hn{l - cy) 
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Fig. 1. — Comparison between different EoS's. Th, n, and Cs, vs for RC (red-long dashed), 
TM (blue-short dashed), ID (green and cyan-dotted), and RP (black-solid). 



-22- 




Fig. 2. — Relativistic shock tube with parallel component of velocity only (PI) with RC 
(red), TM (blue), and ID (green and cyan). 
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Fig. 3. — Relativistic shock tube with parallel component of velocity only (P2) with RC 
(red), TM (blue), and ID (green and cyan). 
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Fig. 6. — Relativistic shock tube with transverse component of velocity (T3) with RC (red), 
TM (blue), and ID (green and cyan). 



